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Abstract 

We investigate both analytically and numerically the renormalization group equations 
in 2D Z{N) vector models. The position of the critical points of the two phase transitions 
for > 4 is established and the critical index v is computed. For N=7, 17 the critical 
points are located by Monte Carlo simulations and some of the corresponding critical 
indices are determined. The behavior of the helicity modulus is studied for N=5, 7, 17. 
Using these and other available Monte Carlo data we discuss the scaling of the critical 
points with and some other open theoretical problems. 
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1 Introduction 



The Berezinskii-Kosterlitz-Thouless (BKT) phase transition was originally discovered in the 
two-dimensional {2D) XY model in the first half of the seventies [H [2]. Since then it was 
realized that this type of phase transition takes place in a number of other models including 
discrete 2D Z{N) models for large enough N and even 3D gauge models at finite temperature 
(see [3111] for a recent study of the deconfinement transition in 3D U{1) lattice gauge theory). 
Here we are interested in the phase structure of 2D Z{N) vector models 0. On a 2D lattice 
A = with linear extension L and periodic boundary conditions, the partition function of the 
model can be written as 



The BKT transition is of infinite order and is characterized by the essential singularity, i.e. 
the exponential divergence of the correlation length. The low-temperature or BKT phase is 
a massless phase with a power-law decay of the two-point correlation function governed by a 
critical index r]. The Z{N) spin model in the Villain formulation has been studied analytically in 
Refs. [6l [7]. It was shown that the model has at least two BKT-like phase transitions when > 
5. The critical index t] has been estimated both from the renormalization group (RG) approach 
of the Kosterlitz-Thouless type and from the weak-coupling series for the susceptibility. It turns 
out that ?7(/3c^'') = 1/4 at the transition point from the strong coupling (high-temperature) 
phase to the massless phase, i.e. the behavior is similar to that of the XY model. At the 
transition point /3c from the massless phase to the ordered low-temperature phase, one has 
ry(/3f ) = A/N\ A rigorous proof that the BKT phase transition does take place, and so 
that the massless phase exists, has been constructed in Ref. [8] for both Villain and standard 
formulations. Monte Carlo simulations of the standard version with A^=6, 8, 12 were performed 
in Ref. [H]. Results for the critical index rj agree well with the analytical predictions obtained 
from the Villain formulation of the model. In Refs. [lOKTT] we have started a detailed numerical 
investigation of the BKT transition in 2D Z{N) models for N = 5, which is the lowest number 
where this transition can occur. Our findings support the scenario of two BKT transitions with 
conventional critical indices. 

Despite of this progress, some theoretical issues remain unresolved. 

• Status of the critical index u. This index governs the behavior of the correlation length 
^, namely it is expected that in the vicinity of the phase transition one has log^ ~ 
O'/iPc ~ P)^ ■ In all numerical studies one postulates the XY value v = 1/2. This value 
appears compatible with numerical data. However, no theoretical arguments exist which 

^See Ref. ^ for a shorter presentation of the resuhs of this work. 
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support this value. Moreover, the study of the Hamiltonian version of Z{N) models via 
the strong coupling expansion combined with the Pade approximant points to a different 
value for A^=5, 6. Only for = 12 one finds z/ 1/2 [6]. 

• The scaling of critical points /3c ' with is still unknown. What is seen from the 
available numerical data is that pi^'' approaches the XY critical point rather fast, probably 

(2) 

exponentially with A^, while /3c scales like some power of A^. 

Most theoretical knowledge about the critical behavior of the XY model comes from the 
solution of the RG equations [2]. Explicit solutions for the RG equations of the arbitrary 
Z{N) model are unknown. Moreover, to the best of our knowledge, even approximate and/or 
numerical solutions have not been obtained. 

One of the goals of the present paper is to fill this gap by constructing both numerical 
and approximate analytical solutions to the RG equations. Another goal is to continue the 
numerical study of 2D Z{N) vector models and to compare results with analytical predictions. 

The paper is organized as follows. In the next section we investigate the system of RG 
equations describing the critical behavior of Z{N) models in the Villain formulation. This 
allows us to calculate the position of critical points for all A^ and to cornpute the critical 
index u. In the section 3 we study numerically the models for N=7 and 17 q We locate the 
transition points and compute some critical indices. Then we present and discuss our results 
for the helicity modulus. Section 4 is devoted to the analysis of the dependence of the critical 
points on A^ and the comparison with the RG prediction. Finally, we summarize our results in 
section 5. 



2 Analysis of RG equations 

The system of RG equations describing the critical behavior of 2D Z{N) models is given by [6] 

where the parameters x, y, z are initially defined as 

X = 7r/3 , ?/ = 27rexp[ — tt^/?] , 2; = 27rexp 



^We chose N ^ 7 since this value interpolates = 6 and iV = 8 which were considered in Ref. ^ and could 
be used as check-points; there is no special reason for the choice of TV = 17, except that we made a balance 
between the need of a value considerably larger than the largest one considered so far in the literature {N = 12) 
and a not too large one to slow down numerical simulations. 
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Figure 1: (Color online) Renormalization group trajectories for = 7 in the vicinity of the 
first phase transition (left) and in the vicinity of the second phase transition (right). The upper 
line shows initial values. 



and a = exp(r) is the lattice spacing [^1. 

This system of RG equations has been obtained for the Villain formulation of the model. 
All available numerical data show that the standard and the Villain formulations belong to the 
same universality class in the case of the XY model. This probably holds also for Z{N) models. 
The only case where universality has been questioned is Z{5), where it was found numerically 
that the helicity modulus has a different behavior in the vicinity of /3c^^ while it jumps to zero 
crossing the critical point in the Villain formulation, it goes to zero smoothly in the standard 
formulation [12]. 

The system ([2]) has two sets of the fixed points. The first one 

x = 2, y = 0, z = 0, (4) 
corresponds to a transition from the massive to the massless phase. The second one 

X = — , y = , z = , (5) 

■^We remind that there are two types of vortices in the Coulomb gas representation of the partition function 
of the 2D Z{N) model residing on the direct and the dual lattices. Variables y and z can be related to the 
corresponding self-energies of these configurations. 
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Table 1: Values of /3c obtained from analytical and numerical solutions of Eqs. ([2]) for various 
values of N. 



N 


6^'^ 

yc an 


5^'^ 

yc num 


yc an 


yc num 


5 


0.726050 


0.741654 


0.841969 


0.853845 


6 


0.739670 


0.747749 


1.212435 


1.219515 


7 


0.740256 


0.747851 


1.650259 


1.659667 


8 


0.740266 


0.747853 


2.155441 


2.167726 


9 


0.740266 


0.747853 


2.727980 


2.743528 


oo 


0.7403 


0.7479 







describes the transition from the massless to the ordered phase. One can see that in the vicinity 
of the first fixed point the behavior of the solution in the r — i- oo limit is the same as in the 
XY model, i.e for /3 ^ /3c^\ x tends to x^o ^ 2, y tends to 0, and for /3 < /3c^\ y tends to 
infinity {z tends to zero in both cases). Thus, /3c^^ may be determined as the value of /3 at 
which the solution goes to the point x = 2, y = z = 0. Clearly, /3c can be determined in a 
similar fashion, studying the system around the second fixed point (3 = ^ (the behavior of y 
and z is reversed). Then, one can easily compute the RG trajectories starting from the initial 
values ([3]). As an example, RG trajectories, projected onto x-y and x-z planes, are shown in 
Fig. [T]for N = 7. Dashed lines show the trajectory at the critical point. Numerical values of 
critical points for various values of N are given in Table [TJ 

To establish how the critical points scale with N and to calculate the critical index z/, it is 
desirable to have an analytical solution for the system ([2]). In the absence of the exact solution, 
we attempt here to construct an approximate solution, valid for large values of A^. Making the 
following change of variables 



1 



X 



y 



exp 



Ar2 ■ 
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Here 



e(r/) 



exp 



^--r+T6 ^ 



and then linearizing the system ([2]) around the first fixed point, one finds 

/ N du ^ dr] ^ 

dr 2 cLT dr 



(6) 



(7) 



is treated as a small function for large A^. This is valid in two cases: 1) A^ ^ 4 and 2) r — t- oo 
if AT > 4. 



5 



To study the system around the second fixed point x = we make the following change 
of variables: 



iV2 



X 



exp 



(9) 



The linearization of the system (|2]) around the second fixed point leads to the same system of 
equations (^^. This means that the solution of the system in terms of ^, u and e will be the 
same as for the first critical point - we only have different initial conditions. 
To the first order in e{r]), one obtains 



X = 2 + - (^0+6: 



^0 



eorfr + Cs . (10) 



Here is the solution of the system in the limit — )■ oo, 

dr + 2 ^° - ^ ' 



^+eo6 = 9{r) 



111 



(12) 



(13) 



(14) 



To find critical points, consider the zero-order approximation (IT^ . In the limit C — )■ (critical 
line of the model) and taking C3 = — 2 In sinh Ci(C2 — tq), the solution can be written as 



^1 describes the first nontrivial corrections and obeys the following linear equation: 

A^2 /■ 

The solution for ,^0 reads 

^0 = 2CiCothCi(r-ro) , 
while the solution for ^1 vanishing for N ^ 00 can be presented as 

= exp(-r7o) / g{T)exp{r]o) dr . 
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We fix To from the requirement that ^o(O) and Co'o(O) satisfy the initial conditions for x and y 
given in ([3]) (this effectively means that ^0 and uq describe the critical trajectory of the XY 
model). Hence, 



, W 1 

27rexp ( — TT 12 

2 V To 



1 

'2^ ' 



To 



2 - nf3^y ■ 



(16) 



From the last equation one finds ~ 0.7403. It is now straightforward to calculate first 
order corrections to this solution. After a long algebra one finds the following equation: 



327r2 exp [-n^P^^'^] = \ - — (/^^'^ - /3f ^) + l-a'^N'' exp 



'0 



4A 



(1) 



(17) 



An approximate solution is given by 



TT 



iV2 



(7r/3,^^ - 2) (2 - 27r + /^^^^tt^) 



exp 



XY 



The same strategy applied for the second fixed point leads to the equation for /3, 



(2). 
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where tq is a solution of the equation 
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These results are summarized in the Table [T] /3c\r? is the approximate value computed 
from (fT7|) (first critical point) and from (fT9|) -( l20|) (second critical point); /3c^num is the numerical 
value computed from the system ([2]). The last row gives the critical values for the XY model: 
the analytical solution /3c = 0.7403 of the linearized system and the numerical solution /3c = 
0.7479 of the exact equations ([2]). This RG prediction should be compared with the Monte 
Carlo result /3c ~ 0.751 for the Villain model p!3]. 

Finally, we would like to gain some information on the critical index v. To compute v one 
has to construct the solution in the region ?/ — > oo for r — > cxd (for the first transition) and 
2; — )■ oo for r — )■ oo (for the second transition). These are the two regions with non- vanishing 
mass gap. The corresponding solutions can be easily derived from the solutions described above. 
It is clear, however, that by the very virtue of that construction, the leading singularity will 
be the same as in the XY model. This means v = Xjl for all large enough X . Therefore, it 
is much more informative if we consider models with not too large and construct fits for r 
from the numerical solutions of the RG equations ([2]) in the given regions. We remind that 



r = loga ~ log^ , 



(21) 



where ^ is the correlation length. Hence, we fit r with the following function 



Table 2: Index v derived from ([2]) and (1221) for the first phase transition. 





A 


B 


V 


/3c 




points 


max (/3c — /3) 


5 


-1.78073 


0.76176 


0.516378 


0.741654 


2.5613 


150 


0.01 


5 


-2.77757 


0.834452 


0.507893 


0.741654 


0.00627118 


100 


0.005 


6 


-3.03443 


0.914498 


0.506457 


0.747749 


0.00744256 


100 


0.01 


7 


-3.09949 


0.928961 


0.504642 


0.747851 


0.00454228 


100 


0.01 


8 


-3.01409 


0.910821 


0.507149 


0.747853 


0.00435793 


100 


0.01 


9 


-3.01271 


0.910531 


0.507189 


0.747853 


0.00422122 


100 


0.01 



Table 3: Index v derived from ([2]) and fl22|) for the second phase transition. 





A 


B 


V 


/3c 




points 


max (/3 — /3c) 


5 


-2.56099 


0.874616 


0.510268 


0.853845 


2.12609 


150 


0.01 


5 


-2.65193 


0.888885 


0.508654 


0.853845 


0.0243011 


100 


0.005 


6 


-3.07805 


1.18864 


0.504577 


1.21951 


0.0212314 


100 


0.01 


7 


-3.04924 


1.38005 


0.505627 


1.65967 


0.00454228 


100 


0.01 


8 


-3.22569 


1.61538 


0.502838 


2.16773 


0.0271673 


100 


0.01 


9 


-3.25341 


1.82077 


0.502766 


2.74353 


0.0212578 


100 


0.01 



In Tables [2] and [3] we give results for N=5, 6, 7, 8, 9. The number of fitting points together with 
the maximal value (/3c — /3) are also shown in the Tables. Fitting points have been distributed 
uniformly between /3c — /3 = and max (/3c — /3). Only central values for all coefficients in both 
Tables are shown. 

Deviations from the central values are in general very small, except for the coefficient A 
for = 5 and when max (/3c — /3) = 0.01. The situation is much improved if one takes 
max (/3c — /3) = 0.005. This can indicate that the scaling region in Z{5) model is somewhat 
narrower than for > 5. We have also tried to fit r assuming the power-like singularity for the 
correlation length. The quality of fits is poor in this case and u acquires a very large value which 
varies with varying A^. This leaves no doubts that the correlation length diverges exponentially 
with u = 1/2 for all A^ > 4 in the vicinity of both phase transitions. 
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3 Numerical results 



We simulated the model defined by Eq. ([T]) using the same cluster Monte Carlo algorithm 
adopted in the Refs. [TOl [11] for the case = 5. We used several different observables to 
probe the two expected phase transitions. In order to detect the first transition {i.e. the one 
from the disordered to the massless phase) we used the absolute value \Ml\ of the complex 
magnetization, 



Ml 

and the helicity modulus [HI [15 



Z^E^^P(^§^^) =1^^!^'^ (23) 



T = {e)-L'/3{s'), (24) 

where e = "^^ij^^ cos {9i — Oj) and s = ^2 X]<ij>^ sin(6'j — 9j) and the notation < ij >x 
means nearest-neighbors spins in the x-direction. 

For the second transition {i.e. the one from the massless to the ordered phase) we adopted 
the real part of the "rotated" magnetization, 

Mr= I Mi I cos W) , 

and the order parameter 



m. 



^ = cos{NiIj) 



introduced in Ref. [16], where tp is the phase of the complex magnetization defined in Eq. ( [2311 . 
In this work, both for N = 7 and = 17, we collected typically 100k measurements for 
each value of the coupling /3, with 10 updating sweeps between each configuration. To ensure 
thermalization we discarded for each run the first 10k configurations. The jackknife method 
over bins at different blocking levels was used for the data analysis. 

In Fig. [2] we show the behavior of the susceptibility xi*^^ = -^^((I^lP) — (I^lI)^) of the 
absolute value of the complex magnetization, which exhibits, for each volume considered, a 
clear peak signalling the first phase transition. The position of the peak in the thermodynamic 
limit defines the first critical coupling, pi^K Fig. [3] shows instead the behavior of versus /3 

(2) 

on various lattice sizes; here the second critical coupling /3c is identified by the crossing point 
(in the thermodynamic limit) of the curves formed by the data on different lattice sizes. 

To determine the first critical coupling l3c^\ we could extrapolate to infinite volume the 
pseudo-critical couplings given by the position of the peaks of xi*^^- However, since the approach 
to the thermodynamic limit is rather slow (powers of logL), we adopted a different method, 
based on the use of the ''reduced fourth-order" Binder cumulant 

^ " 3(|Miiy ^^^^ 
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Figure 2: (Color online) Susceptibility x}^ versus {3 in Z{7) (left) and Z{17) (right) on lattices 
with several values of L. 



the cumulant b[^^^'^ defined as 

(M,j _ {\Mn-{Mn)\') 

-(|M«-(M«)|y' ^''^ 

and the helicity modulus T. We estimated /3c^^ by looking for (i) the crossing point of the 
curves, obtained on different volumes, giving a Binder cumulant versus /3 and (ii) the optimal 
overlap of the same curves after plotting them versus (/3 — /3c) (log L)^/'^, with u fixed at 1/2. 
The method (ii) has been applied also to the helicity modulus T. Our best values for /3c^'' are 

N = 7: = 1.1113(13) , 

N = 17: = 1.11375(250) . 

Then, we performed the finite size scaling (FSS) analysis of the magnetization \Ml\ and the 
susceptibility Xl!^^ at /3c^^ using the following laws: 

IMM'^) = AL-^/\ (27) 

H/3«) = (28) 

where 7/1/ = 2 — and r] is the magnetic critical index. Results are summarized in Tables HI 
O [6] and [71 We observe that the hyperscaling relation ■y/u + 2f3/h' = d = 2 is nicely satisfied 
within statistical errors in both models. 

We can cross-check our determination of the critical exponent rj by an independent method, 
which does not rely on the prior knowledge of the critical coupling, but is based on the con- 
struction of a suitable universal quantity [T71 [11]. The idea is to plot x^^^^^L^"^ versus b[^^^^ 
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Figure 3: (Color online) Behavior of with /3 in Z{7) (left) and Z{17) (right) on lattices with 
several values of L. 

Table 4: Results of the fit to the data of |Mi|(/3c^'') with the scaling law ( 1271) on lattices 
with L ^ Lmin, for N = 7. 



N = 7 





A 




xVd.o.f. 


32 


1.00653(48) 


0.12210(08) 


5.5 


64 


1.00858(70) 


0.12243(12) 


3.7 


128 


1.01074(94) 


0.12277(15) 


2.0 


256 


1.0146(16) 


0.12336(26) 


0.40 


384 


1.0162(22) 


0.12359(34) 


0.19 


512 


1.0177(38) 


0.12381(56) 


0.16 


640 


1.0185(57) 


0.12393(84) 


0.28 



and to look for the value of r] which optimizes the overlap of curves from different volumes. We 
found that, both in Z{7) and Z{17), r] = 1/4 is this optimal value, since it gives the best overlap 
of these curves in the region of values corresponding to the first phase transition, i.e. the lower 
branch of the curves of Fig. |H This result for t] agrees with the determinations rj = 2 — 7/z/ 
from the FSS analysis. 

As for the second critical coupling /3c^\ we used the same method adopted for I3c^\ but 
applied now to and m^. Our best estimates are 

N = 7: /3f ) = 1.8775(75) , = 17 : = 10.13(12) . 

The standard FSS analysis applied to the susceptibility Xl!^^^ of the rotated magnetization Mr 
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Table 5: Results of the fit to the data of \Ml\{I3c ) with the scaling law ( 1271) on lattices 
with L ^ Lrain, for N = 17. 



N =17 




A 




xVd.o.f. 


32 


1.00388(51) 


0.12111(09) 


7.98 


64 


1.00620(69) 


0.12149(12) 


3.58 


128 


1.0089(11) 


0.12191(18) 


1.54 


256 


1.0107(15) 


0.12219(24) 


0.74 


384 


1.0113(24) 


0.12228(36) 


1.36 



Table 6: Results of the fit to the data of (/3c ) with the scaling law (!28|) on lattices 
with L ^ Lmin, for N = 7. 



N = 7 


-^min 


B 


7/z/ 


xVd.o.f. 


32 


0.00558(07) 


1.7402(23) 


2.38 


64 


0.00540(09) 


1.7457(29) 


1.32 


128 


0.00522(12) 


1.7508(38) 


0.68 


256 


0.00518(20) 


1.7520(61) 


0.84 


384 


0.00546(33) 


1.7443(93) 


0.70 


512 


0.00489(52) 


1.760(16) 


0.28 


640 


0.00444(76) 


1.775(25) 


0.0066 



at /SF^ leads to the result for the critical indices 7/z/ given in Tables E and |9]B 

Also in this case the critical index rj can be determined by an independent method, irrespec- 
tively of the knowledge of MrL^/"^ is plotted versus and the value of rj is searched for, 
which optimizes the overlap of data points coming from different volumes. The results we found 
for 7] in Z{7) and Z{17) are in perfect agreement with the theoretical prediction 77*^^^ = 4/A^^ 
(see Fig. [5]). 

Finally, in Fig. Owe present the behavior with /3 of the helicity modulus (1241) . This quantity 
is constructed in such a way that it should exhibit a discontinuous jump (in the thermodynamic 
limit) at the critical temperature separating the disordered phase from the massless one, if the 
transition is of infinite order (BKT). Since the Kosterlitz-Thouless RG equations for the XY 
model [Il[T8l[T9] lead to the prediction that the helicity modulus T jumps from the value 2/(7r/3) 
to zero at the critical temperature, one can check if the same occurs for vector Potts models. In 

'''We do not report in this work the determinations of ji/v by the FSS analysis of the rotated magnetization 
Mr, since they are affected by large statistical and systematic uncertainties. 
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Table 7: Results of the fit to the data of Xl (/^c ) with the scaling law (l28l) on lattices 
with L ^ Lmin, for N = 17. 





iV 


= 17 




-^min 


B 


7/z/ 


xVd.o.f. 


32 


0.00559(08) 


1.7372(26) 


2.7 


64 


0.00532(11) 


1.7453(35) 


0.39 


128 


0.00521(15) 


1.7484(46) 


0.16 


256 


0.00514(23) 


1.7507(71) 


0.15 


384 


0.00522(31) 


1.7483(92) 


0.13 



Table 8: Results of the fit to the data of Xl (Pc ) with the scaling law ( 128|) on lattices 
with L ^ Ljnin, for N = 7. 



N = 7 


Lmin 


A 


7/1/ 


xVd.o.f. 


32 


0.8767(37) 


1.92340(71) 


2.02 


64 


0.8833(47) 


1.92219(87) 


1.41 


128 


0.8858(57) 


1.9217(11) 


1.57 


256 


0.8997(93) 


1.9193(16) 


1.02 


384 


0.916(15) 


1.9166(25) 


0.68 


512 


0.921(24) 


1.9158(39) 


0.98 


640 


0.942(34) 


1.9124(54) 


1.07 



Fig. [6] we plot a red line, representing the function 2/(7r/3); the crossing between this line and 
the curves formed by data points of T approaches indeed when the lattice size increases. 

The knowledge of the behavior with /3 of the helicity modulus provides us with another 
method for the determination of the critical index 77. As shown in Ref. [191 120] (see also 
Ref. [21]), the following relation holds, 

which allows us to get the value of r] at any fixed /3 in the BKT phase, if the value of T at that 
/3 is known. A simple inspection of the behavior of T with /3, shown in Fig. [6], tells us that, for a 
given A^, r] in the BKT phase decreases monotonically from a value compatible with 77*^^^ = 1/4, 
taken at the first critical coupling to a value compatible with r/*^^) = A/N"^, taken at the 
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0.0 - 
0.8 



Z(7) I 



T|=0.25 



• L=32 
■ L=64 

L=128 

' L=256 

- L=384 

• L=512 

> L=1024 



1.0 



1.2 



1.4 



1.6 



2.0 




Figure 4: (Color online) Correlation between x^l ^'^'^L^ ^ the Binder cumulant B)^^^''^ for 
T] = 0.25 in 2'(7) (left) and Z{17) (right) on lattices with L ranging from 128 to 1024. 



AMr) 



Table 9: Results of the fit to the data of Xl Wc ) with the scaling law (!28l) on lattices 
with L ^ Lmin, for = 17^ 



N = 17 


Lfnin 


B 


7/z/ 


xVd.o.f. 


32 


0.9319(47) 


1.98933(89) 


1.65 


64 


0.9408(68) 


1.9878(12) 


1.38 


128 


0.9533(89) 


1.9857(16) 


0.67 


256 


0.954(16) 


1.9856(27) 


0.83 


384 


0.950(24) 


1.9861(40) 


1.09 


512 


0.931(38) 


1.9892(62) 


1.44 


640 


0.911(59) 


1.9925(98) 


2.67 



second critical coupling. In some sense, the iV^ drop of the value of r] at the second transition 
is related to increasing distance between fif'' and 

We studied also the specific heat at the two transitions, finding that, in contrast to the 
case of first- and second-order phase transitions, it does not refiect any nonanalytical critical 
properties at the critical temperatures, thus confirming that only BKT transitions are at work 
here. 
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Figure 5: (Color online) Correlation between MrLP^'^ and for r] = 4/7^ in Z{1) (left) and 
Tj = 4/17^ in Z{n) (right), on lattices with L ranging from 128 to 1024. 




Figure 6: (Color online) Helicity modulus versus /3 in Z{5), Z{7) and Z{17) on lattices with 
various sizes. 



4 Behavior with N of the critical couphngs 

The results of this work and those available in the literature allow us to make some considera- 
tions about the behavior with of the critical coupling Pc^^ and Examining Table [Done 
concludes that formulae ffTK]) and f lT^ give the correct qualitative prediction for the scaling of 
the /3c ' with in the Villain formulation. Namely, 

• Pc^^ converges to XY value very fast, like exp(— aA^^) 

• /Sc^-* diverges like A^^ . 

One should expect that the standard vector Z{N) model possesses similar scaling, probably up 
to 0{N) corrections. One could try therefore to fit available Monte Carlo data for /3c ' with 
formulae f|T8|) and f|T9|) modified to account for such corrections. 
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Table 10: Summary of the known values of the critical couplings Pc^^ and for 2D Z{N) 
vector models. 







PP 


Reference 


5 


1.0510(10) 


1.1048(10) 


m 


6 


1.11012(74) 


1.4257(22) 


[9J 


7 


1.1113(13) 


1.8775(75) 


this work 


8 


1.11907(88) 


2.3480(22) 


[9] 


12 


1.11894(88) 


5.0556(128) 




17 


1.11375(250) 


10.13(12) 


this work 


oo 


1.1199(1) 


oo 


m 



In Table [TO] and in Figs. [71 |8]we summarize the present knowledge about the position of the 
critical points for 2D Z{N) vector models. 

We can see from Fig. [7] that the approach of Pc^^ to the N oo limit, corresponding to the 
2D XY model, is indeed very fast. Introducing corrections into the scaling formula ( 1T8|) . we 
may conjecture the following general behavior: 

/3^^) = A-{BN^ + CN + D)exp 

Indeed, attempts to fit data with this formula give values for A and E compatible with 
However, as is seen from Table [IHl the combination of our Monte Carlo determinations and 
those of Ref. [9] leads to a fake non-monotonic approach to in marked contradiction with 
our RG analysis. Looking separately at the two data sets given in Table [TOl one can see that 
each set indeed satisfies monotonicity (taking into account the error bars in the case of data from 
Ref. [9]). We think, therefore, that the non-monotonicity of the combined set can be explained 
with the different systematics affecting the determinations of pi^^ by the two collaborations. 
On our side, we have determined the location of Pc^^ by two different methods, using several 
observables and working on larger lattices, up to L = 1024, thus making us very confident on 
the reliability of our results. In conclusion, a reliable check of the formula f[30l) . based on all 
data given in Table [TOl is not possible and that formula remains a conjecture. 

In the case of /3c , instead, one can verify that, for example, the following extension of f[T^ 

/3f = ^ + BN + C + De-'T (31) 

fits the data rather well, although with a high x^/d.o.f. probably reflecting the different sys- 
tematics mentioned above, 

A = 25.89(115) , B = -0.040(23) , C = 0.38(9) , D = 227.7(634) , = 12-69 . 



(30) 
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Figure 7: (Color online) Behavior with of the known values of the critical coupling Pc^'' in 
2D Z{N) vector models. The horizontal dashed line represents the critical coupling of 2D XY, 
/3xy=l-1199, taken from Ref. |22j. 



The value of A is close to Svr, as is expected. 

In general, however, we have to conclude that even if the scaling formulae (130|) and (1311) 
might be correct, more data for different N are needed to determine all coefficients. 

5 Summary 

In this paper we have studied the RG equations describing the critical behavior of 2D Z{N) 
vector models in the Villain formulation. The main original results are 

• the RG trajectories in the vicinity of both phase transitions, 

• the critical points as functions of A^, 

• the index u, which turns out to be equal to 1/2 for all > 4. 

The numerical part of the work has been devoted to verify these theoretical expectations: 
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Figure 8: (Color online) Behavior with of the known values of the critical coupling in 
2D Z{N) vector models. 



• We have determined numerically the two critical couplings of the 2D Z{N = 7, 17) vector 
models and given estimates of the critical indices rj at both transitions. Our findings 
support for all > 5 the standard scenario of three phases: a disordered phase at 
high temperatures, a massless or BKT one at intermediate temperatures and an ordered 
phase, occurring at lower and lower temperatures as increases. This matches perfectly 
with the iV — )• oo limit, i.e. the 2D XY model, where the ordered phase is absent or, 
equivalently, appears at /3 — > oo. 

• We have found that the values of the critical index rj at the two transitions are compatible 
with the theoretical expectations. 

• The index v also appears to be compatible with the value 1/2, in agreement with RG 
predictions. 

On the basis of this study and taking into account previous works [9l [11] , we are prompted 
to conclude that 2D Z{N) vector models in the standard formulation undergo two phase tran- 
sitions of the BKT type. Furthermore, the standard and the Villain formulations belong to the 
same universality class with v = 1/2, rj^^^ = 1/4 and r/*-^^ = A/N'^. 
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Considering the determinations of the critical couphngs as a function of A^, we have cal- 
culated the leading dependences using RG equations and conjectured the approximate scaling 
for /3c ' ( A^) in the standard version. The existing numerical values and their accuracy are not 
sufficient, however to reliably check the conjectured formulae. 
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